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We study a recently introduced model of one-component glass-forming liquids whose constituents 
interact with anisotropic potential. This system is interesting per-se and as a model of liquids 
like glycerol (interacting via hydrogen bonds) which are excellent glass formers. We work out the 
statistical mechanics of this system, encoding the liquid and glass disorder using appropriate quasi- 
particles (36 of them). The theory provides a full explanation of the glass transition phenomenology, 
including the identification of a diverging length scale and a relation between the structural changes 
and the diverging relaxation times. 



The study of associated liquids like glycerol as glass 
formers has a long and rich history [l[ , but until now the 
role of the anisotropic hydrogen bonds, while clearly im- 
portant in frustrating crystallization, has not been made 
explicit. Recently a model of one component liquids with 
anisotropic interaction potential was introduced Q, to- 
gether with numerical simulations in two-dimensions that 
demonstrated clearly the importance of the anisotropic 
interaction in frustrating crystallization and allowing the 
formation of a glassy state of matter. This model is 
important in stressing the fact that even simple one- 
component liquids may not crystallize if the local symme- 
try of the interaction potential does not perfectly match 
the symmetry of the equilibrium crystal. It is worthwhile 
therefore to analyze further this example of glass forma- 
tion and put it in the general context of the glass tran- 
sition. In this Letter we present a theory of this model, 
constructing its statistical mechanics and providing an 
understanding of the phenomenology of its glass transi- 
tion, including an identification of a diverging length and 
explaining the diverging time scales. Our analysis allows 
putting this interesting example of glass formation on the 
same footing as other classical glass formers such as bi- 
nary mixtures with centralpotentials 0, 0] , stressing the 
generality of the approach [J 0| and of the glass transition 
phenomenon at the same time. 

Particles of mass m in this model interact via 
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where Vij is the distance between the two particles i and 
j. The first term on the RHS of ^ is the standard 
isotropic Lennard- Jones potential 
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whereas the anisotropic part of the potential is given by 
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h{x) = {l-x^f for \x\ < 1 ;h{x) =Ofor|.T| > 1 . 



Here 9i {0j) is the included angle between the relative 
vector Tij = ri — rj and a unit vector ut (uj) (referred 
to below as 'spin') which represents the orientation of the 
axis of particle i {j). The function h{{9 — 9q)/9c) (with 
6'o = 126° and 9c = 53.1°) has a maximum at 6* = 6*0, 
and thus is a favored value of 9i. Thus the anisotropic 
term in the potential favors structures of five-fold sym- 
metry. The parameter A controls the tendency of five- 
fold symmetry, and therefore of the frustration against 
crystallization. The units of mass, length, time and tem- 
perature are m, a, t = a^Jm/e and e/kB, respectively, 
with kB being Boltzmann's constant. 

According to the numerical simulations presented in 
0, for A < 0.6 this system crystallizes upon reducing the 
temperature. The ground state crystal has an elongated 
hexagonal structure with anti-ferromagnetic ordering of 
the spins Ui , but the actual crystal that is obtained upon 
cooling is a "plastic crystal" with hexagonal spatial order 
but with spin disorder. For A > 0.6 the system fails to 
crystallize upon cooling. The relaxation times were mea- 
sured by monitoring the rotational autocorrelation func- 
tions Cnit) = {l/N)J2Mt{t) ■ u,{0)) which was fitted 
to a stretched exponential form C/?(t) cx cxp[—{t/Ta)^]. 
For A = 0.6 the relaxation is of Arrhenius form with a 
constant value /3 w 0.95 for T > r,„ = 0.46, but /3 was 
fit separately for every temperature T < where it 
decreases with temperature. The relaxation times were 
fitted to a Vogel-Fulcher law t„ = tq exp[DTo/(r - Tq)] 
which involves fitting the three free parameters tq = 0.61, 
D = 7.4 and Tq = 0.099 (in addition to p). We repeated 
the simulations of this model using Monte-Carlo meth- 
ods in N-P-T ensemble finding results in agreement 
with the MD simulations of in the same ensemble. 

To construct the statistical mechanics of this system 
we recognize that the potential energy between any pair 
of particles depends on their spin orientations. In Fig. 
[l]wc present the three potentials between two particles, 
depending on the orientation of their spins relative to the 
inter-particle vector distance: lowest in energy (in blue 
continuous line) is the case for which both have a favored 
spin orientation; middle, in green dashed line (high, in 
red dotted line) is the potential when one (none) of the 
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FIG. 1: (Color online). Potential curves for particle-pairs 
with two spins, one spin or no spin in favored position (blue 
continuous line, green dashed or red dotted line respectively). 
Inset: the measured energies of particle pairs, falling in three 
distinct ranges with gaps between them, allowing us to define 
the quasi-particles. The peak in each colored curve corre- 
sponds to the minimum in the main figure. 



spins are in a favored orientation. One sees that the 
minima of these potentials occur with significant gaps in 
their energies, allowing us to now measure the average 
energy of pairs of particles as a function of temperature. 
These averages fall in three distinct ranges, such that the 
range of variation of each energy is much smaller than the 
gaps between the energies, see inset in Fig. [TJ This al- 
lows us to proceed to define quasi-species. We denote the 
three effective energies below as 2Eb 2Eg and 2Er respec- 
tively. Note that the spin orientations involved in each 
such mean energy can fluctuate within a temperature- 
dependent range of angles. For the temperature range of 
interest the range of angles is in a sector of about 60°, 
but as this range of angles determines the degeneracies 
that enter the statistical mechanics below, we need to 
reconsider it carefully as explained in the sequel. Next 
consider an n-star, which by definition is a given particle 
decorated by the n inter-particle vector distances (edges) 
to its n neighbors, see for example Fig. [2] Each such edge 
is colored according to the spin orientations. We denote 
by fc the number of red, green and blue edges such 
that n = i -\- j + k. It turns out that in the temperature 
range of interest (0 < T < 0.5), in an overwhelming ma- 
jority of n-stars (more than 98%) the central particle has 
a spin orientation that is favorable with respect to two 
of its edges (this is of course the maximal value, which 
is favored by energy considerations). Therefore we take 
a-priori j + k > 2^ neglecting the very small number of 
instances where this docs not hold. The energy of such 




FIG. 2: (Color online). An example of an n-star with n — 
5,1 — 2, j — 2 and k — 1. The central particle has a spin 
with favored orientation with respect to edges 1 and 2. Thus 
these edges can be either blue or green, and this central spin 
cannot be favored with respect to any other edge. In the 
interesting range of temperatures we observe 36 n-stars with 
4:<i + j + k <6. 

an n-star (referred to as a quasi-particlc) is computed as 

E^jk = iEr + jEg + kEb , (4) 

where k < 2. Note that since the energies on the RHS 
of Eq. (jl]) depend on temperature, so does the energy 
of the quasi-particles. Notwithstanding, in the inter- 
esting temperature range the temperature dependence 
is weak; we take the energies of the quasi-particles as 
T-independent (we used as half the energy of a parti- 
cle pair Er = -0.2187, Eg = -0.5645, Eb = -1.5105). 
The degeneracy gijk of the energy level (number of quasi- 
particles with the same energy) is 

Since the central particle always has a spin in favorable 
orientation to two edges, each of these two edges must be 
blue or green; the first factorial is the number of possible 
choices of blue edges; one way if there are two (or none) 
of them and two ways if there is one of them. Once these 
two edges are determined, there remain i + j + k — 2 
edges to choose the i red from, giving rise to the second 
factorial. This completes the degeneracy due to color. 
Next we count the number of spin orientations. There 
are 2 — k green edges due to the central particle and i red 
edges, giving us i — k+2 unfavorable spin orientations and 
j + 2A: — 2 favorable ones. The number of ways to orient 
the unfavorable spins is 4»-'=+2 and the number of ways to 
orient the favorable spins is 2^+2*^-2 ^ rpj^^ number 4 stems 
from the fact that the two favored orientations occupy 
an angular sector of 2 x 60°, leaving us with 4 sectors of 
60° for the unfavorable orientations. The fact that the 
central particle can emanate at most two favored edges 
means that we have a constraint X^ijfcO + 2fc)cy7; < 4, 
where Cijk is the mol-fraction of quasi-particles having i. 
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FIG. 3: (Color online). Direct numerical simulation of the 
constraint 

j and k edges of the right color. In practice, as mentioned 
above, the inequality can be swapped with an equality 

+ 2fc)cyfc = 4 . (6) 

ijk 

To satisfy this constraint exactly we need to consider the 
temperature dependence of the spin-fluctuation sectors, 
since when these change, so do the assignments of (ijk). 
In Fig. [3] we show the LHS of Eq. [6] for two fixed spin- 
fluctuation sectors (upper, red squares, 69°, lower, green 
diamonds, 50°, and middle, blue circles, variable sector 
of width w = (68. ST + 41.2)°). In aU cases the sum 
was measured via Monte Carlo simulations in the range 
0.05 < T < 0.5. The quahty of the constraint using the 
variable spin-fluctuation sector is obvious. The decrease 
at high temperatures is due to the increased fluctuations 
in the spin orientations and in the energies of the quasi- 
particles, inducing changes in the degeneracies and in 
the {ijk) assignments. We thus use this temperature 
dependent width of the spin-fluctuations to assign the 
quasi-particle index (ijk) in all our simulations. 
Now write the partition function of the system: 

z(r,A(r)) = ^5y,e-'^^-- 6-^(^+2'=) . (7) 

ijk 

The Lagrange multiplier A is introduced to insure that 
the constraint ([6]) is satisfied. In terms of the partition 
function the mol-fraction of quasi-particles is 

Z{T,X{T)) ■ ^' 

Substituting Eq. dH) in Eq. §^ we compute A(T) for 
each temperature, and then compute the mol-fraction 
Cijk- For presentation and comparison with numerical 
simulations it is advantageous to bunch groups of Cijk 
together. One bunching is in the three groups obtained 




FIG. 4: (Color online). Comparison of the direct numerical 
simulations to the theoretical prediction. Shown are the mol- 
fractions Ck with fc = (dashed-dotted red line), k — 1 (green 
dashed line) and k = 2 (continuous blue line). One sees the 
point of departure of the direct numerical simulations from 
equilibrium (the point of jamming) which is estimated to be 
about T = 0.17. Inset: same concentrations on logarithmic 
scale. 

with k ~ 0,1,2. In Fig. [Jwe present a comparison of the 
theory to the simulation for the mol-fractions of quasi- 
particles with fc = 0,1,2. We note that the agreement 
is excellent down to T « 0.17 where the simulation gets 
jammed. This observation is in agreement with Q who 
estimated the glass transition temperature to be about 
0.15 on the basis of the "divergence" of relaxation times. 
Wc note that the statistical mechanics predicts the pre- 
cise spin statistics of the glassy jammed state, since the 
mol-fractions freeze at a "Active temperature" T « 0.17 
that can be read directly from Fig. 21 We do not need 
to measure relaxation times to see where the system falls 
out of equilibrium; it is obvious directly from Fig. 21 

In addition, the statistical mechanics predicts two 
"transitions" when the mol-fractions of quasi-species 
with fc = and then with fc = 1 become small. The 
glass transition (jamming) occurs visibly when the mol- 
fraction Ck=o becomes small. Wc read a second transi- 
tion when Ck=i becomes exponentially small. This im- 
plies that the ground state consists solely of fc = 2 quasi- 
particles, in perfect agreement with the existence of the 
a crystalline ground state with anti-ferromagnetic order. 
Thus the offered statistical mechanics explains very well 
the phenomenology of this system. Note that both these 
"transitions" refer to finite systems, where sufficiently 
small concentrations mean effectively zero concentration. 

The greatest challenge for the statistical mechanics 
is whether it can also predict the measured relaxation 
times. Jamming is caused by the rapid reduction in the 
mol-fraction of some spin configurations, leading to a loss 
of ergodicity. To see this consider the quasi-particles with 
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FIG. 5: (Color online). Comparison of the measured spin 
rotational relazxation times Ta (points) with the relaxation 
time predicted by Eq. (|10|) (continues line). 

k = 0. These are the highest in energy and accordingly 
their mol-fraction goes to zero first when the tempera- 
ture cools down. Using then the mol-fraction Ck=o hi 
comparison with the area per particle a = A/N, we form 
a length scale ^ according to 

^ V^/ck^oiT) . (9) 

When Ck=o — *■ the length scale ^ — *■ oo, defining regions of 
increasing size that are jammed and therefore contribut- 
ing to an increasing relaxation time. Spin relaxations are 
dominated by correlated stringy (1-dimensional) chains, 
and we estimate the number of quasi-particlcs involved, 
N*, as N*{T) w ^{T). The relaxation time is determined 
by the free energy barrier, and denoting by A* the chemi- 
cal potential per quasi-particle we write [5j, |9| 
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Note that this prediction differs essentially from the 
Adam-Gibbs formula Q in the sense that it does not pre- 
dict a divergent Tq at any finite temperature, but rather 
an enormous increase in Tq, when Cg=o ~^ exponentially 
fast. The statistical mechanics does not allow Cg=o = 0. 



Of course, in any finite system "exponentially small" mol- 
fractions can be actually zero, and the relaxation time 
can be effectively infinite. The theory does not recognize 
however a sharp transition in the thermodynamic limit. 

The comparison of the prediction ()10p to the measured 
values of the relaxation times is shown in Fig. [5l We 
note both the excellent agreement and the fact that tq 
is of the order of unity, as expected in the limit T ^ oo 
where the relaxation time should be the particle vibra- 
tion time. We again draw the attention of the reader 
to Fig. 3 in where a fit to the relaxation time is 
achieved, using similar ideas. In the problem there the 
relaxation was configurational rather than via an internal 
variable as here, and typically relaxation events spanned 
2-dimensional correlated domains, making the free en- 
ergy barrier proportional to The equal usefulness of 
the ideas used, with the only change in the identification 
of the quasi-species, their degeneracy and the constraints 
on the statistical mechanics gives us hope that the ap- 
proach is quite general and can be applied to glass form- 
ing systems of very different nature. Whether or not such 
a computer-assisted statistical mechanics can be applied 
to 3-dimensional glass formers is a question that must 
await future research. 
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